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Abstract 


Geometric  approaches  for  filling-in  surface  holes  are  introduced  and  studied  in  this 
paper.  The  basic  idea  is  to  represent  the  surface  of  interest  in  implicit  form,  and 
fill-in  the  holes  with  a  scalar,  or  systems  of,  geometric  partial  differential  equations, 
often  derived  from  optimization  principles.  These  equations  include  a  system  for 
the  joint  interpolation  of  scalar  and  vector  fields,  a  Laplacian-based  minimization,  a 
mean  curvature  diffusion  flow,  and  an  absolutely  minimizing  Lipschitz  extension.  The 
theoretical  and  computational  framework,  as  well  as  examples  with  synthetic  and  real 
data,  are  presented  in  this  paper. 

Keywords:  Inpainting,  variational  formulations,  interpolation,  surface  holes,  scalar  and 
vector  fields,  Laplacian,  mean  curvature,  absolute  minimizing  Lipschitz  extension. 
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1  Introduction 


Inpainting  is  a  term  used  in  art  to  denote  the  modification  of  images  (painting,  pho¬ 
tographs,  etc)  in  a  form  that  can  not  be  detected  by  an  ordinary  observer.  It  normally 
refers  to  the  filling-in  of  regions  of  missing  information  or  the  replacement  of  regions  by 
a  different  kind  of  information.  This  is  a  very  important  topic  in  image  processing,  with 
applications  including  image  coding  and  wireless  image  transmission  (e.g.,  recovering  lost 
blocks),  special  effects  (e.g.,  removal  of  objects),  and  image  restoration  (e.g.,  scratch  re¬ 
moval).  The  basic  idea  behind  the  computer  algorithms  that  have  been  proposed  in  the 
literature  is  to  fill-in  these  regions  with  available  information  from  their  surroundings. 
This  information  can  be  automatically  detected  as  in  [12,  25],  or  hinted  by  the  user  as  in 
more  classical  texture  filling  techniques  [22,  26,  38].  Several  names  have  been  used  for  this 
filling-in  operation,  including  disocclusion  in  [8,  32],  or  inpainting  in  [11,  12,  13].  In  the 
context  of  this  paper,  and  following  [12] ,  we  shall  refer  to  it  as  digital  inpainting. 

It  turns  out  that  images  are  not  the  only  kind  of  data  where  there  is  a  need  for  digital 
inpainting.  Surfaces  obtained  from  range  scanners  often  have  holes,  regions  where  the  3D 
model  is  incomplete.  The  main  cause  of  holes  are  occlusions,  but  these  can  also  be  due  to 
low  reflectance,  constraints  in  the  scanner  placement,  or  simply  lack  of  sufficient  coverage 
of  the  object  by  the  scanner.  This  is  frequently  observed  in  the  scanning  of  art  pieces  [31], 
and  is  in  part  due  to  the  fact  that  complicated  geometry  has  a  lot  of  self-occlusions  and 
details.  Art  pieces  also  impose  significant  restrictions  on  the  scanner  placement.  Holes  are 
also  observed  in  common  scenarios  where  LADAR  data  is  collected  (e.g.,  a  house  behind 
an  occluding  tree),  and  in  all  the  major  areas  where  range  scanners  are  used.  With  the 
increasing  popularity  of  range  scanners  as  3D  shape  acquisition  devices,  with  applications 
in  geoscience,  art  (e.g.,  archival),  medicine  (e.g.,  prohestetics),  manufacturing  (from  cars 
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to  clothes),  and  defense  (e.g.,  LADAR),  it  is  very  important  to  be  able  to  inpaint  this 
missing  information.  This  is  often  needed  for  post-processing  as  well  as  for  presentation. 
It  is  the  goal  of  this  paper  to  present  a  framework  for  inpainting  these  surface  holes. 

Our  work  is  inspired  by  the  one  reported  in  [21],  and  it  is  presented  as  an  alternative  to 
this  method.  This  pioneering  work  addressed  the  problem  of  hole  filling  via  isotropic  diffu¬ 
sion  of  volumetric  data  (that  is,  iterative  Gaussian  convolution  of  some  distance  function 
to  the  known  data).  The  approach  proposed  by  these  authors  addresses  holes  with  com¬ 
plicated  topology,  a  task  very  difficult  with  mesh  representations.  The  reader  is  directed 
to  this  paper  for  an  excellent  and  detailed  description  of  the  nature  of  holes  in  scanning 
statues  and  for  a  literature  review  in  the  subject.  We  should  only  note  that  most  algo¬ 
rithms  on  reconstructing  surfaces  from  range  data  are  point-cloud  reconstruction  based 
and  treat  holes  as  regions  with  low  sampling  density,  thereby  interpolating  across  them 
[2,  6,  10,  18,  24,  27].  Of  course,  these  algorithms  often  do  not  distinguish  between  a  real 
hole  in  the  data  and  one  due  to  the  lack  of  sampling,  and  equally  fill  or  fail  to  fill  both 
cases  in  the  same  fashion.  Other  point-cloud  methods  evolve  a  surface  over  time  until  it 
approximates  the  data  [17,  42,  44],  or  fit  a  set  of  3D  radial  basis  functions  to  the  data, 
compute  a  weighted  sum  of  them  and  use  a  level  set  of  this  last  function  as  reconstructed 
surface  [23,  14].  Mesh  based  methods  for  surface  reconstruction  [39,  20,  41]  can  perform 
hole  filling  as  a  post- process  or  integrate  hole  filling  into  surface  reconstruction  [20].  One 
of  our  proposed  models  is  closely  related  to  the  one  presented  in  [19]  (and  of  course  to 
our  previously  introduced  3D  surface  inpainting  model  [40]),  where  the  authors  use  the 
Willmore  flow  with  a  finite  element  implementation.  In  contrast  with  their  work,  our 
model  works  on  implicit  surfaces,  thereby  allowing  for  more  complicated  hole  topologies, 
and  also  naturally  leads  to  systems  of  low  order  differential  equations. 
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The  first  algorithm  here  proposed  is  an  extension  of  our  previous  work  on  image  in¬ 
painting  [7,  8,  12]  (see  also  [11,  13,  16,  32,  34,  37]).  In  particular,  we  show  how  to  adapt 
the  variational  formulation  we  presented  in  [7,  8]  to  the  problem  of  surface  hole  filling.  As 
in  [21],  the  use  of  volumetric  data  (that  is,  the  surface  is  represented  as  the  zero  level-set 
of  a  function)  brings  us  topological  freedom.  In  contrast  with  [21],  we  use  a  system  of 
coupled  anisotropic  (geometric)  partial  differential  equations  designed  to  smoothly  con¬ 
tinue  the  isophotes  of  the  embedding  function,  and  therefore  the  surface  of  interest  (as 
the  zero  level  isophote).  These  equations  are  based  on  the  geometric  characteristics  of 
the  known  surface  (e.g.,  the  curvatures),  and  as  [21],  are  applied  only  at  the  holes  and  a 
neighborhood  of  them  (being  these  equation  anisotropic  and  geometry  based,  they  lead 
to  a  slightly  slower  algorithm  than  the  one  reported  in  [21],  as  expected  with  geometric 
flows).  A  preliminary  version  of  this  (first)  model  was  presented  in  [40].  We  formalize 
this  and  improve  it  here  with  an  automatic  initialization  method.  This  initialization  is 
based  on  the  computation  of  a  conical  neighborhood  T  of  the  known  part  of  the  surface, 
call  it  S,  where  the  distance  function  is  uniquely  attained.  Thereby  we  can  define  the 
signed  distance  function  ds  and  then  Vds  is  the  extension  of  the  unit  normal  to  S  to  a 
neighborhood  of  it.  This  construction  also  helps  us  to  label  both  parts  of  the  surface  as 
interior  and  exterior,  and  this  is  useful  in  this  first  method. 

We  also  develop  additional  curvature  based  hole  surface  inpainting  methods.  The 
first  of  them  is  based  on  a  variational  model  which  integrates  the  Laplacian  of  a  distance 
function  (i.e.,  a  function  which  satisfies  |VD|  =  1,  and  D  =  ds  in  the  conical  neighborhood 
T),  in  a  open  set  containing  the  hole.  Recall  that  the  Laplacian  of  the  distance  function 
gives  the  mean  curvature  of  its  level  sets.  The  second  method  is  more  heuristic  and  is 
based  on  the  diffusion  of  a  function  which  represents  mean  curvature  of  level  sets  of  an 
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underlying  implicit  function. 


Finally,  we  also  present  simpler  methods  based  on  the  Laplace  equation  and  the  so- 
called  AMLE  model,  which  permit  to  reconstruct  a  function  which  is  distance-like  near  the 
known  part  of  the  surface  and  whose  zero  level  set  can  be  interpreted  as  the  reconstructed 
surface.  If  our  interest  is  just  to  find  a  smooth  reconstruction,  this  approach  may  be 
sufficient.  If  one  wants  a  reconstruction  which  is  based  on  minimizing  mean  curvature,  it 
can  serve  as  an  initialization. 

These  algorithms,  except  the  one  based  on  curvature  diffusion  which  is  less  reliable, 
exhibit  a  similar  behavior  in  reconstructing  surface  holes  for  synthetic  and  real  data.  As 
mentioned  above,  the  reconstructions  based  on  the  Laplace  or  AMLE  equation  can  be  used 
as  initializations  for  the  curvature  based  approaches.  Describing  and  studying  all  these 
techniques  provides  a  comprehensive  understanding  of  the  different  possible  frameworks 
for  filling-in  surface  holes. 

The  remainder  of  this  paper  is  organized  as  follows.  Section  2  describes  our  variational 
approach  for  image  inpainting  in  any  dimension  N  and  its  adaptation  to  the  reconstruction 
of  surface  holes.  In  Section  3  we  present  the  two  additional  curvature  based  approaches: 
a  variational  one  minimizing  the  integral  of  the  absolute  value  (or  a  power  of  it)  of  the 
mean  curvature,  and  an  heuristic  one,  based  on  the  diffusion  of  curvature.  Section  4 
describes  two  simple  methods  for  surface  reconstruction  based  on  Laplace  equation  and 
the  so-called  AMLE  model.  Section  5  describes  the  numerical  algorithms  used  in  our 
computations.  In  Section  6  we  present  some  numerical  experiments  on  hole  filling  obtained 
with  the  algorithms  previously  introduced.  Finally,  in  Section  7,  we  summarize  the  main 
conclusions  of  the  paper. 
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2  Joint  interpolation  of  vector  fields  and  gray  levels  and  its 


application  to  surface  inpainting 

Let  us  describe  the  variational  approach  to  filling-in  by  joint  interpolation  of  vector  fields 
and  gray  values  which  was  introduced  in  [7,  8] .  Our  purpose  is  to  adapt  it  to  the  problem 
of  hole  filling  on  surfaces. 

Let  Q  be  a  hyper-rectangle  in  MN  and  P  an  open  bounded  subset  of  Q  with  smooth 
boundary.  Suppose  that  we  are  given  an  image  uo  :  Q  \  P  — * >  M,  where  P  denotes  the 
closure  of  P.  Using  the  information  of  uq  on  Q  \  P  we  want  to  reconstruct  the  image  uq 
inside  the  hole  of  missing  information  P.  In  our  context,  the  function  uq  is  an  implicit 
representation  of  the  known  data.  In  [7,  8]  we  proposed  to  fill-in  the  hole  P  (on  images) 
using  both  the  gray  level  u  and  the  vector  field  of  normals  9  to  the  level  sets  of  the  image 
outside  the  hole.  This  permitted  to  design  energy  functionals  which  minimize  a  power  of 
(mean)  curvature  and  to  write  them  in  terms  of  the  pair  of  variables  ( u,9 ).  This  is  the 
approach  that  we  shall  explore  next  with  the  purpose  of  interpolating  holes  in  surfaces. 

We  denote  by  LP(Q),  1  <  p  <  oo,  the  space  of  (measurable)  functions  /  :  Q  — »  M  such 
that  [q  \f(x)\p  dx  <  oo.  By  L°°(Q)  we  denote  the  space  of  bounded  functions  /  :  Q  — >  JR. 

Let  P  be  an  open  subset  of  Q  with  smooth  boundary  such  that  P  CC  P.  The  band 
around  P,  used  to  fill  it  in,  is  the  set  B  =  P  \  P.  To  fill-in  the  hole  P  we  use  the 
information  of  uq  contained  in  B,  mainly  the  gray  level  uq  and  the  vector  field  of  normals 
(i.e.,  the  gradient  directions)  to  the  level  sets  of  uq  in  B,  which  we  denote  by  9q.  We 
assume  that  Oq  is  a  vector  field  with  values  in  MN  satisfying  Oq(x)  ■  Vuo(x)  =  |Vrto(x)| 
and  |#o(z0|  ^  1-  Intuitively,  this  amounts  to  say  that  Qq(x)  =  when  Vuo(x)  /  0, 

and  we  give  freedom  to  6q(x)  when  Vuq(x)  =  0  requiring  only  that  |0q(^)|  <  1-  The  basic 
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goal  then  is  to  extend  in  a  smooth  way  the  pair  ( uq,6q )  from  the  band  B  =  17  \  17  to  a 
pair  of  functions  (it,  0)  inside  17.  For  that  we  attempt  to  continue  the  isosurfaces  of  uq 
(i.e.,  the  hypersurfaces  [«o  =  A]  :=  {x  :  uo(x)  =  A}  or,  more  generally,  the  boundaries  of 
the  level  sets  [uq  >  A] ,  A  £  1R )  in  B  inside  17  by  taking  into  account  the  principle  of  good 
(interpreted  here  as  smooth)  continuation.  The  energy  functional  proposed  in  [7,  8]  was 
based  on  the  following  principles: 

a)  We  constrain  the  solution  (it,  8)  to  coincide  with  the  data  on  the  band  B.  For  that  we 
ask  that  u  =  uq  on  B.  The  vector  field  8  should  also  satisfy  |0|  <  1  on  17  and  should  be 
related  to  u  by  8  •  Vu  =  |Vu|,  i.e.,  we  impose  that  8  is  the  vector  field  of  directions  of  the 
gradient  of  u.  This  implies  that  8{x)  =  6q{x)  on  the  points  x  €  B  where  Vuo(x)  /  0  and 
we  leave  freedom  to  6{x)  on  the  points  where  x  £  B  where  X7uq(x)  =  0  besides  satisfying 
that  |0(x)|  <  1. 

b)  We  impose  that  the  vector  field  8q  in  the  band  B  is  smoothly  continued  by  8  inside  17. 
Note  that  if  8  are  the  directions  of  the  normals  to  the  level  hypersurfaces  of  u,  then  div(0) 
(a  possible  measure  of  smoothness  of  the  vector  field)  is  the  mean  curvature.  The  smooth 
continuation  of  the  levels  sets  of  uq  inside  17  is  imposed  by  requiring  that  div0  £  Lp(!7), 
i.e.,  J~  |div  8 |p  dx  <  oo.  For  consistency  we  shall  require  that  8q  is  such  that  div  8q  £  LP(B) 
and  assume  that  (ito5^o)  maY  be  extended  to  a  a  pair  of  functions  (u,0)  in  the  class  of 
functions  where  we  minimize  our  energy  (see  (1)). 


Based  on  these  basic  principles,  we  proposed  in  [7,  8]  to  interpolate  the  pair  ( 8 ,  u)  in 


17  by  minimizing  the  functional 


Minimize  f  |div(#)|p(7  +  /3|Vk  *  u|)dx 

Jn 


\9\  <  1,  II  u  ||oo<  M, 


|  Vit|  —  9  •  V'u  =  0  in  17, 


(1) 


u  =  uo  in  B,  9  •  vn  =  9q  ■  in  <917, 

where  p  >  1,  7  >  0,  /?  >  0,  k  denotes  a  regularizing  kernel  of  class  Cl  such  that  k(x )  >  0 
a.e.,  and  M  =  ||uo||l°°(_b)  :=  supxe^  |uo(x)|.  As  usual,  denotes  the  outer  unit  normal 
to  <917.  Let  us  note  that  if  u  is  the  characteristic  function  of  a  set  A  C  ]R3  (i.e.,  u(x)  = 
Xa{%)  =  1  if  x  G  A  and  =  0  otherwise)  with  smooth  boundary  and  9  is  a  smooth  extension 
of  the  unit  normal  to  dA,  then  f~  |div(0)|p|Vu|dx  =  fgA  |H|pdS,  where  H(x)  is  the  mean 
curvature  of  dA  and  dS  denotes  the  surface  area.  The  convolution  of  Vu  with  the  kernel 
k  is  done  for  technical  reasons,  it  permits  to  prove  the  existence  of  a  minimum  for  (1) 
[7,  8],  though  we  can  dismiss  it  in  practice.  Finally,  let  us  also  note  that  the  constant  7 
has  to  be  >  0,  it  implies  an  Lp  bound  on  div  9,  and  this  is  useful  to  prove  that  the  limits 
of  minimizing  sequences  satisfy  the  constraint  |  Vw|  —  9  ■  Vu  =  0  in  17  [7,  8].  We  refer  to  [8] 
for  a  detailed  theoretical  analysis  of  this  formulation  and  its  approximation  by  smoother 
functionals.  Let  us  finally  note  that  the  use  of  smooth  continuation  principles  based  on 
powers  of  mean  curvature  as  smoothness  measure  was  proposed  in  [33]  with  the  purpose  of 
image  segmentation  with  depth  (reconstructing  then  the  occluded  boundaries),  and  used 
as  a  base  for  variational  approaches  of  image  disocclusion  in  [7,  8,  32], 
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2.1  Surface  inpainting 


Let  us  now  describe  how  to  adapt  the  above  formulation  to  inpaint  (fill-in)  holes  (or  gaps) 
on  surfaces  S,  which  we  assume  to  be  embedded  in  JR3.  To  avoid  any  confusion  with  our 
previous  use  of  the  word  hole,  let  us  use  the  word  gap  of  the  surface.  Assume,  to  fix  ideas, 
that  S  is  a  smooth  compact  connected  surface,  and  At  is  a  part  of  S  which  is  unknown  or 
could  not  be  obtained  during  scanning  (or  is  damaged  and  needs  to  be  replaced).  Let  us 
identify  S  with  its  known  part.  Let  us  choose  a  bounding  box  Q  in  JR3  strictly  containing 
the  surface  gap  M.  and  part  of  S  (see,  for  instance,  Figure  8).  Let  dJv[  be  the  boundary 
of  the  gap.  Even  if  J\4  is  unknown,  its  relative  boundary  in  S  is  known.  Let  A"  be  a 
neighborhood  of  S  n  Q  such  that 

F  =  {x  G  Q  :  d(x,  S  n  Q)  <  ad(x,  dM)},  0  <  a  <  1. 

We  assume  that  F\ (S  n  Q)  consists  of  two  connected  components,  which  can  be  identified 
as  the  two  sides  of  the  surface  S  (see  Figure  1).  The  information  derived  from  the  region 
F  is  considered  reliable  and  we  impose  it  as  a  constraint  in  our  reconstruction.  Let  df  (x) 
be  the  distance  of  a  point  x  £  T  to  S  n  Q.  By  changing  the  sign  of  dT  in  one  of  the  sides 
of  the  surface  we  may  define  the  signed  distance  function  to  S  n  Q  in  T  (take  it  positive 
inside  and  negative  outside).  We  denote  it  by  df(x),  or  simply,  by  ds.  The  vector  field  in 
F- 

N(x)  =  Vds(x), 

is  an  extension  of  the  unit  normal  vector  field  on  S  n  Q  to  its  neighborhood  F .  Again,  we 
consider  this  information  as  reliable  and  it  will  be  used  as  a  constraint. 

To  adapt  functional  (1)  to  surface  hole  reconstruction  we  must  make  explicit  the  hole 
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(a)  (b) 

Figure  1:  a)  Section  of  the  surface  S  with  the  hole  M  and  the  neighborhood  T.  b)  Sign 
assignment  to  the  two  faces  of  S.  The  hole  17  is  defined  as  the  ball  minus  T.  The  band  B 
is  defined  as  Q  \  17. 

17  and  the  functions  (no,  Oo)  which  are  the  known  data  on  a  neighborhood  of  17.  To  define 
the  hole  we  take  a  ball  Q  (or  any  open  set  homeomorphic  to  a  ball)  such  that  Q  CC  Q  and 
containing  the  boundary  of  the  gap  dA4  in  its  interior.  We  define  the  hole  17  by  removing 
from  Q  the  points  of  T .  We  take  the  band  B  =  Q  \  17  (notice  that  Q  is  identified  to  the 
set  17  of  Section  2). 


Figure  2:  If  the  ball  is  not  large  enough,  we  may  have  difficulties  in  the  assignment  of 
labels.  Notice  that  if  we  extend  the  ball,  the  labelling  can  be  done  in  a  consistent  way. 
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We  then  consider  uq  :  Q  \  If  — >  1R  a  characteristic  function,  that  is,  a  binary  function 
taking  values  0  and  1.  The  values  uq(x)  =  1  and  uq(x)  =  0  represent  the  points  which  are 
interior,  respectively,  exterior,  to  S.  Recall  that  we  are  assuming  that  (<Sn  Q)  consists 
of  two  connected  components  which  represent  the  two  sides  of  the  surfaces.  We  then  label 
these  two  sides  with  the  values  uq  =  1,  representing  the  inner  part  of  the  surface,  and 
uq  =  0,  representing  its  outer  side.  The  set  [uq  =  1]  HJ7  (respectively,  the  set  [uq  =  0]  H^7) 
is  represented  in  Figure  l.b  with  the  sign  +  (resp.,  the  sign  -).  By  propagation,  we  extend 
this  labelling  to  the  rest  of  Q  \  If,  knowing  it  already  in  T.  For  simplicity,  we  assume 
that,  by  eventually  extending  the  ball,  this  labelling  can  be  done  in  a  consistent  way:  we 
assume  that  Q  can  be  chosen  so  that  dQ  is  divided  into  two  labelled  regions  consistent 
with  the  labels  in  B  \  (S  n  Q)  (see  figures  l.b  and  2).  We  call  A  the  set  of  points  x  in 
Q\Q  such  that  uo(x)  =  1,  hence  u o(x)  =  xa(x)-  In  this  case,  by  minimizing  (1),  we  want 
to  reconstruct  the  set  A  inside  the  hole  f I  knowing  the  set  outside  If. 

Recall  that  we  minimize  (1)  by  solving  the  gradient  descent  equations  (9),  (10)  pre¬ 
sented  below,  using  the  numerical  approach  described  in  Section  5.1,  and  we  need  the 
initial  conditions  for  u  and  9  in  Q  =  11  U  B,  so  that  8  ■  Vu  =  |Vit|.  For  that,  we  define 
Uinit  in  Q  =  fl  U  B  as  the  extension  of  uq  inside  11  by  a  geodesic  propagation,  so  that  u;nit 
takes  values  0  and  1.  Then  we  define  Ds  as  the  signed  distance  to  <9[uin;t  =  0]  (negative 
in  ['Uinit  =  0]  and  positive  in  =  1])  so  that,  by  construction,  Ds  is  an  extension  of  ds 
(recall  that  ds  is  defined  in  T).  We  take  the  vector  field  0jn;t  =  V D s  in  Q  =  OU B  and  we 
define  9q  =  0jnjt  in  B.  Observe  that  Vrto  =  vSsnQ-;  where  v  denotes  the  inner  unit  normal 
to  S  and  dsnQ  is  the  Hausdorff  measure  on  S  n  Q.  We  have  9m •  Vujnit  =  |  V rtinit  |  - 
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3  Alternative  curvature  based  approaches 


Following  the  description  in  the  introduction,  we  now  present  alternative  filling-in  ap¬ 
proaches  for  the  problem  of  filling  holes  in  three  dimensions. 

3.1  Energy  in  terms  of  distance  functions 

Recall  that  if  S'  is  a  smooth  manifold  of  class  C 2,  then  the  signed  distance  function  D  to 
S'  is  also  of  class  C 2  in  a  neighborhood  of  S' .  The  vector  field  \7D  is  an  extension  of  the 
unit  normal  to  S'  and  satisfies  |VU|  =  1.  The  operator  A D{x)  =  div  VD(x)  represents  the 
sum  of  the  principal  curvatures  of  the  isosurface  [D  =  D(x)\  :=  {y  G  Q  :  D(y)  =  D(x)}. 
When  we  look  at  this  as  a  function  in  JR3,  the  distance  function  is  Lipschitz  and  it  satisfies 
\VD\  =  1  in  the  viscosity  sense.  The  isosurfaces  may  develop  singularities  and  the  only 
thing  we  can  expect  is  that  the  mean  curvature  is  a  Radon  measure.  Indeed,  recall  that 
the  mean  curvature  of  a  polyhedral  surface  is  a  Dirac’s  measure  concentrated  at  the  edges, 
and  the  signed  distance  function  may  have  such  type  of  singularities.  We  shall  assume  that 
the  signed  distance  D  to  the  surface  S  is  such  that  AD  G  JA(Q )  where  M(Q)  denotes 
the  space  of  Radon  measures  in  Q  [1],  We  define 

W{Q)  =  {«£  L\Q)  :  Vu  g  L\Q),  A  u  G  M(Q)}. 

We  propose  to  fill-in  the  three  dimensional  holes  via  the  minimization  of  the  functional 

Min{D  G  W(Q),  \VD\  =  l,D  =  dsmF}fQ \AD(x)\dx-  (2) 

This  energy  integrates  the  mean  curvature  on  the  isosurfaces  of  D.  Due  to  the  singularities 
of  the  isosurfaces  of  D,  the  integral  of  a  power  of  the  mean  curvature  with  an  exponent 
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p  >  1  may  be  infinite.  Let  us  observe  that  problem  (2)  has  a  minimizer  as  soon  as 
the  admissible  set  is  nonempty,  and  we  assume  that  this  is  the  case.  Details  on  the 
implementation  of  the  proposed  approach  are  presented  in  Section  5.2. 

3.2  Curvature  diffusion  and  distance  reconstruction 

We  also  present  studies  based  on  diffusion  of  the  mean  curvature  (see  also  [43]  for  related 
work  based  on  the  linear  Poisson  equation).  For  convenience,  let  us  write  Qjr  ■=  Q\J7. 
We  propose  to  diffuse  the  mean  curvature  of  S  and  then  reconstruct  the  surface  with  the 
prescribed  curvature,  that  is,  we  propose  to  solve  the  system  of  PDEs 


ut  =  An; 


in  [0,  oo)  x(Q\f) 


ut  =  \Vu\  (div(^)-u;)  in  [0,  oo)  x  (Q  \  F) 
with  the  following  boundary  conditions  on  dJ-  n  Q: 


(3) 


u j  =  A ds  on  8T  l~l  Q, 


Vu  ■  v  =  Vds  ■  v  on  dT  n  Q, 


(4) 


where  ds  denotes  the  signed  distance  to  S  and  v  denotes  the  outer  unit  normal  to  dJ- . 
Observe  that  we  did  not  write  the  Dirichlet  boundary  condition  u  =  ds  because  it  is  not 
possible,  in  general,  to  impose  it  in  a  classical  sense  to  the  equation  for  u  (mean  curvature 
type  equation).  Let  us  comment  on  the  boundary  conditions  used  on  dQ?  \  dJ7.  First  of 
all,  we  observe  that  the  ideal  scenario  would  be  to  consider  Q  =  1RN  and  solve  the  system 
of  PDEs  (3)  in  [0,  oo)  x  (MN  \JT)  with  the  boundary  conditions  (4),  but  this  is  impossible 
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at  the  numerical  level.  For  that  we  modify  the  boundary  conditions  on  dQ?  \  dT ’,  we  do 
linear  extrapolation  of  u  and  u  along  the  normal  to  the  level  sets  of  u.  Details  on  the 
implementation  of  the  proposed  flow  are  presented  in  Section  5.2. 

4  The  Laplace  and  the  Absolute  Minimizing  Lipschitz  ex¬ 
tension  interpolation 

In  [15]  we  studied  and  classified  interpolation  algorithms  which  satisfy  a  reasonable  series  of 
axioms  in  terms  of  the  solution  of  a  partial  differential  equation.  Two  particular  examples 
are  the  Absolutely  Minimizing  Lipschitz  Extension,  denoted  as  AMLE  in  the  sequel,  and 
the  Laplacian  interpolation.  We  now  discuss  the  applicability  of  AMLE  and  the  Laplace 
equation  to  the  problem  of  filling-in  surface  holes. 

The  Laplacian  interpolation  is  based  on  solving  the  PDE 


—A u  =  0  in  Qjf,  (5) 

with  specified  boundary  data  on  dQp.  Indeed,  boundary  data  is  only  known  in  8T  n  Q, 
where  we  should  impose  that  u  =  ds.  Thus,  a  reasonable  assumption  would  be  to  consider 

—  =  0  in  8Qt  \  dT  (6) 

where  v  denotes  the  outer  unit  normal  to  dQ ?  \  dT .  In  some  sense,  from  the  theoretical 
point  of  view,  the  lack  of  knowledge  of  boundary  conditions  for  u  in  dQp  \  dT  is  an 
obstacle  to  the  possibility  of  using  (5)  to  reconstruct  the  surface  S  fl  Q  (which  is  defined 
as  d[u  >0]).  In  spite  of  this,  we  can  solve  (5)  with  the  boundary  condition  (6),  or  use, 
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instead  of  (6),  a  linear  extrapolation  of  u  along  its  gradient  direction  on  dQj?  \  T.  The 
results  obtained  with  this  last  boundary  condition  are  presented  in  Section  6.  We  should 
of  course  mention  that  this  approach  is  closely  related  to  the  work  in  [21],  based  on  linear 
diffusion. 

The  AMLE  interpolation  ([4,  5])  is  based  on  solving  the  PDE 

(D2u{Xu)  ,Xu)  =  0  in  Qr.  ( 7 ) 

with  boundary  data  on  dQ ?  (here  Xu  and  D2u  denote  the  gradient  and  the  Hessian  matrix 
of  u,  respectively,  so  that  in  Cartesian  coordinates,  ( D2u  (Xu)  ,  Xu)  =  Yaj=i  ;/,■  g.r  I!','  g‘,‘  )■ 
As  a  short  review  on  (7),  let  us  recall  that  this  equation  can  be  solved  on  general  domains 
and  boundary  data,  in  particular  the  data  can  be  given  in  a  finite  number  of  surfaces, 
curves  and/or  points.  Indeed,  we  may  assume  that  the  boundary  data  ip  G  Lipg(Qjr) 
where 

Lipd{Qr)  =  \g  £C{dQr)  :\\\g\\\=  sup  <  oo]  , 

[  x,yedQjr  dgQ{X,  y)  J 

and  dQr(x,y)  is  the  geodesic  distance  between  x  and  y  in  Qjr,  i.e.,  the  minimal  length  of 
all  possible  paths  joining  x  and  y  and  contained  in  Q ?  [29] . 

Let  us  recall  that,  if  X  is  an  open  set  or  a  smooth  manifold  in  MN ,  IT1’00 (A)  (resp. 
W1,P{X))  denotes  the  space  of  functions  u  e  L°°{X )  (resp.  u  6  LP(X))  such  that  Vn  G 
L°°(X)  (resp.  Xu  G  IX (X)).  By  W01’oo(A)  (resp.  Wq’p(X)),  we  denote  the  closure  in 
Wl,°°(X)  (resp.  in  Wl'p(X)),  of  the  smooth  functions  with  compact  support  in  X. 

Existence  and  uniqueness  of  viscosity  solutions  for  the  AMLE  model  (7)  with  boundary 
data  p  G  Lipg(Qjr)  was  proved  by  Jensen  [29].  Moreover,  he  proved  that  the  viscosity 
solution  of  (7)  is  an  absolutely  minimizing  Lipschitz  extension  of  <p,  i.e.,  u  G  n 
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C(Qjr)  and  satisfies 


I|Vm||loo(Q/;RJV)  <  ||Vu>||Loo(Q/;RJV)  (8) 

for  all  Q'  C  Qj:  and  w  such  that  u  —  w  G  Wq’°° (Q') .  Let  us  add  that  the  AMLE  model 
was  introduced  by  Aronsson  in  [4,  5]  as  the  Euler-Lagrange  equation  of  the  variational 
problem  (8)  (which  can  be  interpreted  as  the  limit  as  p  — >  oo  of  the  variational  problems 
||Vu||  lp(Q';Rn)  <  1 1  ^ 1 1  lp  (Q'-,RN )  f°r  Q'  Qt  and  w  such  that  u  —  W  G  Wq'P{Q') 
[9,  29]).  The  above  results  were  extended  in  [29]  to  the  case  of  continuous  boundary 
data  ip  €  C(dQj?),  and  Jensen  proved  that  in  that  case,  the  AMLE  is  locally  Lipschitz 
continuous  in  Q^p  [29]. 

Concerning  boundary  conditions,  the  same  remarks  we  made  for  the  Laplace  equation 
(5)  can  be  done  here,  that  is,  boundary  data  are  only  known  in  dJ-  n  Q  where  we  should 
impose  that  u  =  ds  (by  the  results  in  [30],  there  exist  absolutely  minimizing  Lipschitz 
extensions  of  ds\gj?rQ  which  satisfy  (7),  but  there  is  no  uniqueness  result  for  them).  From 
the  theoretical  point  of  view,  the  lack  of  knowledge  of  boundary  conditions  on  dQp  \  dJ-  is 
an  obstacle  to  the  possibility  of  using  (7)  to  reconstruct  the  surface  SnQ  (which  is  defined 
as  d[u  >0]).  To  overcome  this  difficulty,  we  propose  to  use  either  Neumann  boundary 
conditions  (6),  or  to  use  linear  extrapolation  of  the  values  of  u  along  the  direction  Vu  on 
QQt  \  The  experiments  displayed  in  Section  6  were  obtained  with  the  last  method. 
Details  on  the  implementation  of  both  Laplacian  and  AMLE  interpolation  methods  are 
presented  in  Section  5.3. 
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5  Numerical  considerations 


We  now  present  some  basic  concepts  related  to  the  numerical  implementation  of  the  differ¬ 
ent  filling-in  models  described  above.  For  later  use  let  us  define  the  backward  and  forward 
derivative  operators  by: 

=  u(i  +  l,j,k)  -  u{i,j,k)  8~xu{i,j,  k)  =  u(i,j,k )  -  u{i  -  1  ,j,k) 

$£>(*> .7>fe)  =  uih3  +  M)  -  u{i,j,k)  6 ~2u(i,j,k)  =  u(i,j,k )  -  u{i,j  -  1  ,k) 

Si3u{hj,k)  =  u(i,j,k+  1)  —  u(i,j,  k)  8~zu{i,j,k)  =  u(i,j,k )  -u(i,j,A;-l) 

where  indexes  ( i,j,k )  denote  the  voxel  location.  We  use  the  superscript  c  to  denote 

centered  differences  in  space  derivatives,  and  we  define  Vcu  =  (#^14,  <5£2u,  8°  u)  where 
8rxi u(i.  j.  k )  =  |  (u(i  +  1,  j,  k)  —  u{i  —  1,  j,  k )),  with  analogous  expressions  for  8^2u{i,j,  k ), 
8^u(i,j,k). 

5.1  Joint  interpolation  of  vector  fields  and  gray  levels 

To  minimize  the  functional  (1)  we  use  the  steepest  descent  method.  If  we  denote  the 
energy  by  E(6,u),  the  steepest  descent  equations  are 

0t  =  -VeE(6,u)  (9) 

and 

ut  =  -VuE(6,u )  (10) 

in  (0,  oo)  x  0,  supplemented  with  the  corresponding  boundary  data  and  initial  conditions. 
The  constraints  on  (9,  u )  can  be  incorporated  either  by  penalization  or  by  projecting  onto 
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them  after  each  time  step.  We  tested  both  methods  in  an  implicit  (also  in  an  explicit) 
in  time  discretization  of  (9),  (10).  Let  us  explain  in  some  detail  the  implicit  in  time 
implementation  of  (9),  (10)  with  the  constraint  9-Vu  =  |Vu|  incorporated  by  penalization. 
The  constraint  u  =  uq  in  B  is  also  incorporated  by  penalization.  Thus  we  consider 

E(u,9)  =  f  |div(#)|p(7  +  /3|Vk  *  u|)dx  +  ??  f  (|Vu|  —  9  •  Vu)  +  A  [  (u  —  uo)2dx  (11) 

which  corresponds  to  the  energy  (1)  plus  two  penalization  terms  for  the  constraints 
9  ■  Vu  =  |Vit|,  and  u  =  uq  on  B,  with  77,  A  >  0.  At  the  discrete  level,  the  integral 
is  replaced  by  a  sum,  the  gradients  are  discretized  with  forward  derivatives  and  the  di¬ 
vergence  with  backward  derivatives  (or  conversely,  gradients  with  backward  derivatives, 
while  divergences  with  forward  ones).  We  follow  this  structure  since  at  the  continuous 
level  the  divergence  is  the  dual  operator  of  the  gradient,  while  at  the  discrete  level  the 
divergence  discretized  with  backward  (resp.,  forward)  derivatives  is  the  dual  of  the  gradi¬ 
ent  with  forward  (resp.,  backward)  derivatives.  Thus,  we  use  as  discrete  approximations: 
Vu  =  (<5+  u,  8+2u,  (5+ it),  and  similarly  for  Vk*u,  while  div  9  =  Xu=i  8^.9\  if  9  =  {9±,  $2,  83). 
To  simplify  our  notation,  let  us  write  g(9)  :=  /3|div(0)|p,  h(u )  :=  7  +  /3\Vk  *  u\.  Then 

VeE{9,  u)  =  —pV[h(u)  |div(0)|p-2div(0)]  -  77 Vu  =  0  (12) 


and 

VuE(9,u )  =  -div  (k*  (g(6)  ^  ^  ^ -??  div  ^  +77  div  0+2A(u-uo)xb  =  0,  (13) 
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both  equations  holding  in  17.  As  we  mentioned  above,  the  gradients  are  discretized  with 
forward  derivatives,  while  backward  ones  are  used  to  discretize  the  divergences.  To  solve 
equations  (9)  and  (10),  we  use  an  implicit  discretization  in  time.  To  be  precise,  we  write 

VeE{9,6',u,v)  =  —pV[h(u)(e  +  |div(0/)|p_2)div(0)]  -  r/Vu,  (14) 


and 


VuE{0,d',u,v)  =  -div  I  k  *  (g(0) 


Vk  *  u 


\J e  +  |Vk  *  v|2 
+V  div(6*)  +  2A(u  -  u0)xb- 


—  r]  div 


Vu 


\/e  +  |Vv|2 


(15) 


Then,  we  use  the  discretization  in  time  given  by 


9n+ 1  -  =  -A tX7eE(dn+\6n,un,un),  (16) 


and 

un+l  _un=  _A tVuE(9n+1,9n+1,un+\un).  (17) 

Finally,  we  make  the  change  of  variables  ^n+1  :=  Qn+1  —  9n ,  vn+1  :=  un+1  —  un  and  we 
have 


Cn+1  =  -A  tVeE{C+1  +  9n,9n,un,un),  (18) 

vn+i  =  _A  tVuE(9n+1,9n+1,vn+1  +  un,un).  (19) 

In  practice  we  solve  (18)  in  17  with  the  boundary  condition 

0n+i  .  ^  =  9q  .  Qn 
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Now,  since  9n+1  ■  =  &n  '  we  may  write  this  boundary  condition  as 


'an 


£n+1  •  vn  =  0  on  dn. 


The  boundary  condition  for  (19)  is  (see 


k  *  (s(*”+1)  >  +  ->  f  ,V"l+1„,,'l  ■  ri  =  -A '  ri  on  an. 


\fe  +  \Vk  *  un |2 


\/e  +  |  Vun|2 


which  is  written  in  terms  of  vn+l  =  un+1  —  un  as 

5fc(0n+1,un+1,un)  =  T]d0  ■  J1  -  Bk(9n+1,un,un), 

where  =  t .  («(<>"+■) ^ ^  +  n  ( -^-S— )  •  A  Then  we  use 

a  conjugate  gradient  method  to  solve  (18)  and  (19).  The  constraint  \6\  <  1  is  incorporated 
by  renormalizing  6n  (when  \6n\  >  1)  after  each  time  step.  The  constraints  on  HuHoo  can 
be  also  introduced  after  each  time  step.  In  spite  of  the  penalization  term,  the  relationship 
6  •'Vu  =  |Vu|  is  (numerically)  lost  and  we  reinforce  it  after  a  certain  number  of  time  steps. 

We  can  also  set  rj  =  0  and  incorporate  the  constraint  that  |Vu|  =  6  ■  Vu  by  projecting 
onto  it  after  each  time  step.  We  also  tested  this  both  in  a  time  implicit  and  explicit 
discretization  of  equations  (9)  and  (10).  After  each  time  step  of  6  and  u  we  redefine 


)  = 


OjiJJ)  +  aVu(i,j,l) 
max(  1,  \9(i,  j,  l )  +  aVu(i,j,  Z)|) 


for  some  a  >  0.  As  it  has  been  shown  in  [28]  this  is  a  good  way  of  imposing  that  |#|  <  1 
and  9  ■  Vu  =  |Vu|.  We  have  found  quite  similar  results  using  both  described  methods  for 
imposing  the  constraint. 
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In  our  experiments,  we  take  k  a  Gaussian  kernel  with  small  variance,  say  one  or  two 


pixels.  In  practice,  one  can  also  dismiss  the  kernel  k.  The  initial  conditions  for  u  and  8 
are  taken  as  we  explained  at  the  end  of  Section  2.1  so  that  8  -X7u  =  |Vu|  in  Q  :=  11  U  B. 
We  could  also  construct  an  initialization  (u,  8)  using  the  solution  Laplace  equation,  or 
AMLE  (see  Section  4).  As  regards  the  parameters  in  functional  (11),  in  our  experiments 
we  choose:  (3  =  1,  A  =  100,  7  =  1  and  7  =  10  (although  we  have  experimented  with 
7  =  0,7  =  0.5  and  7  =  1  and  the  results  are  always  good).  The  parameter  e  in  equations 
(14)  and  (15)  is  set  to  0.1. 

Finally,  let  us  observe  that,  since  both  constraints  u  =  u$  in  B  and  8  ■  Vu  =  |Vu|  in 
fl  have  been  imposed  by  penalization,  this  gives  some  robustness  to  noise,  but  this  can  be 
also  achieved  with  a  previous  smoothing  of  the  surface  near  the  hole. 


5.2  Curvature-base  approaches 

We  use  the  steepest  descent  method  also  to  minimize  the  functional  (2),  and  we  solve 


A  =  -A 


(20) 


In  order  to  satisfy  the  constraint  |VH|  =  1  we  make  use  of  the  PDE  that  computes  the 
signed  distance  function  [35]: 


A  =  -sign(A(|VA  -  1) 


(21) 
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Then,  as  a  numerical  approach  to  minimize  (2)  we  combine  (20)  and  (21)  at  each  time 
step.  To  go  from  Dn  to  Dn+l  we  first  solve 

/  A  Dn 

D*  =  D"  -  At  A  — - ; - 

\\ADn\  +  e 

using  the  standard  discretization  for  the  Laplacian  (Aw  =  J2i= l  ^x  ^x  w  f°r  anY  function 
w)  and  then  solve 

Dn+i  =  D*  _  Atsign(D*)(|VD*|  -  1) 

using  an  upwind  scheme  for  the  gradient  magnitude  [35,  36].  The  precise  numerical  scheme 
is: 


Dn+i  =  d*  -  At  [max(sign(T)*),  0)(|V+H*|  -  1)  +  min(sign(D*),  0)(|V_D*|  -  1)] 


where 


V+D* 


V~D* 


-  3 

J2  max(^T)*,  O)2  +  m.m(5+D*,  0) 
.1=1 

r  3 


max('5ilD*’°)2  +  min(<5T.H*,  0)2 


(22) 


(23) 


Li=l  J 

As  boundary  conditions,  we  use  Dn  =  D*  =  ds  on  dT  fi  Q  and  linear  extrapolation  of  Dn 
and  D*  along  its  gradient  direction  on  8Qf  \  dT . 

Even  if  not  fully  theoretically  justified,  unless  we  work  in  a  small  neighborhood  of  the 
surface  D  =  0,  a  similar  scheme  may  be  used  to  minimize  the  functional 


Min 


[D:  \VD\  =  1,  A De  L2(Q),  D  =  ds  in  JQ 


\AD(x)\2dx.  (24) 


In  order  to  solve  the  system  (3),  we  solve  first  the  equation  =  A oj  using  an  explicit 
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Euler  scheme  and  the  standard  discretization  for  the  Laplacian.  Then,  we  reconstruct  the 
distance  function  by  solving  ut  =  |Vit|  (div  For  that,  we  use  the  following 

explicit  scheme  ([36]): 


un+1  =  un  —  At  [|V+rtn|  max(w  —  div  V^i/\  0)  +  | V  un\ min(w  —  div  V+Kn,  0)] 


where 


div  V+un  =  (6  5  5  ) 


where  e  >  0  is  a  small  parameter,  and  |V+'Un|,  |V  un\  are  given  as  in  (22),  (23). 


5.3  The  Laplace  equation  and  the  AMLE 

The  Laplace  equation  (5)  is  solved  by  computing  the  steady  state  of  equation  ut  =  A u. 
This  PDE  is  discretized  with  an  implicit  (or  explicit)  Euler  scheme  in  time  and  the  standard 
discretization  for  the  Laplacian.  In  case  of  the  implicit  Euler  scheme,  the  corresponding 
linear  system  is  solved  using  the  conjugate  gradient  method.  As  boundary  conditions  we 
use  u  =  ds  on  dJ-DQ  and  linear  extrapolation  of  u  along  its  gradient  direction  on  dQ^\dT 
(one  can  also  use  the  Neumann  boundary  conditions  (6)).  The  same  boundary  conditions 
are  used  for  the  AMLE  equation  whose  numerical  scheme  we  describe  now. 

The  AMLE  equation  (7)  is  also  solved  computing  the  steady  state  of  its  associated 
evolution  problem.  Following  [15],  we  use  a  nonlinear  overrelaxation  method  (with  an 
implicit  Euler  scheme  in  time  and  centered  differences  for  space  derivatives).  Thus,  we 
compute  un+l  at  ( i,j,k )  by: 


un+l  =un_w 


(l  +  2Af)|WjV  -un\Vcu*\2  -  Ai£ij=i 


u%.Xj6^.u* 


(1  +  2Af)|Vc'U*|2  +  e 


(25) 
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where  the  relaxation  parameter  is  0  <  w  <  2  (we  use,  in  practice,  w  =  1.5). 

When  solving  (25),  voxels  are  visited  in  a  lexicographic  order.  Let  us  denote  by  A 
the  lexicographic  order  for  indexes  in  Z3.  Now,  when  computing  the  value  of  ttn+1(i,  j,  k) 
we  use  u* ,  defined  below,  which  is  a  function  that  takes  the  same  values  as  un  or  un+l 
depending  at  the  pixel  location  ( p,q,r )  where  we  evaluate  it: 

{un(p,q,r )  if  (i,j,k)  A  (p,q,r), 

(26) 

un+l  (p,  q,  r )  otherwise, 

Finally,  we  define  u*.x  as 

8~.5£.u*  —  2 u*  if  i  =  j, 

5cXi6cXju*  otherwise. 

This  iterative  scheme  (25)  gives  the  solution  of  (7). 

6  Experimental  results 

We  now  show  experimental  results  illustrating  the  filling-in  techniques  here  proposed. 

6.1  Simple  geometric  objects 

First,  we  present  experiments  of  geometric  objects  done  with  the  different  methods  dis¬ 
cussed  above:  the  joint  interpolation  of  vector  fields  and  gray  levels  (abbreviated  JIVFGL) , 
minimization  of  the  absolute  value  of  the  Laplacian  of  the  distance  function  (and  also  the 
case  of  power  2),  curvature  diffusion,  the  AMLE,  and  Laplace  equations.  The  images  in 
our  experiments  have  been  rendered  using  the  AMIRA  Visualization  and  Modeling  System 
[3], 
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Figure  3  shows  a  pyramid  with  a  hole  and  its  corresponding  reconstructions;  the  basis 


of  the  pyramid  is  a  square  whose  sides  have  length  61  voxels,  the  height  is  of  30  voxels  and 
the  hole  is  parallel  to  the  basis  going  from  voxels  of  height  8  to  voxels  of  heigth  12,  both 
included.  Figure  3(a)  displays  an  ideal  pyramid  from  which  we  cut  a  piece,  the  pyramid 
with  the  hole  is  displayed  in  Figure  3(b).  Notice  that  the  boundary  of  the  hole  looks 
rounded,  this  is  only  due  to  our  display  and  does  not  correspond  to  the  data.  Figure  3(c) 
shows  its  reconstruction  by  joint  interpolation  of  vector  fields  and  gray  levels  (JIVFGL), 
i.e.,  functional  (1)  with  N  =  3.  Figure  3(d)  shows  the  result  obtained  minimizing  the 
absolute  value  of  the  Laplacian  of  the  distance  function,  i.e.,  (2).  Figure  3(e)  shows  the 
result  obtained  minimizing  the  square  of  the  Laplacian  of  the  distance  function,  i.e.,  (24). 
Figure  3(f)  shows  the  result  obtained  with  curvature  diffusion  (3).  Finally,  Figures  3(g) 
and  3(h)  display  the  results  obtained  solving  the  AMLE  and  Laplace  equations  in  3D, 
respectively.  Note  how  the  reconstructions  obtained  with  the  model  JIVFGL  and  the 
square  of  the  Laplacian  for  example  manage  to  fill-in  a  relatively  large  hole.  The  others 
do  a  decent  work  that  can  certainly  be  used  as  a  very  good  initial  condition  for  the  best 
models  to  refine.  For  a  better  display  of  the  results  we  show  in  Figure  4  an  orthoslice  of 
each  of  the  reconstructions,  corresponding  to  a  plane  parallel  to  the  basis  on  the  middle 
of  the  hole  (voxels  located  at  heigth  10). 

Figure  5  displays  a  torus  with  a  hole  and  its  corresponding  reconstructions;  the  torus 
is  generated  by  the  revolution  of  a  circle  of  radius  r  =  10  voxels  along  a  circle  of  radius 
R  =  20  voxels,  the  hole  extends  over  an  angle  of  c  =  35°  (approximately).  Figure  5(a) 
displays  the  ideal  torus  from  which  we  cut  a  piece.  Figure  5(b)  shows  the  gap  in  the 
torus,  notice  again  that  the  boundary  of  the  hole  looks  rounded,  this  is  only  due  to  our 
display  and  does  not  correspond  to  the  data.  Figure  5(c)  shows  its  reconstruction  by  joint 
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interpolation  of  vector  fields  and  gray  levels  (JIVFGL).  Figure  5(d)  displays  the  result 
obtained  minimizing  the  absolute  value  of  the  Laplacian  of  the  distance  function,  i.e. ,  (2). 
Figure  5(e)  displays  the  result  obtained  minimizing  the  square  of  the  Laplacian  of  the 
distance  function,  i.e.,  (24).  Figures  5(f)  and  5(g)  display  the  results  obtained  solving 
the  AMLE  and  Laplace  equations  in  3D.  Note  that  the  best  results  are  obtained  with  the 
JIVFGL  model  and  with  (24).  For  a  better  display  of  our  results,  we  show  in  Figure  6  an 
orthoslice  of  each  reconstruction  through  a  plane  passing  by  the  middle  of  the  torus. 

6.2  Experiments  with  Michelangelo’s  David 

For  this  real  data,  with  the  purpose  of  adapting  it  to  our  algorithm,  the  data,  originally 
given  as  a  triangulated  surface,  was  converted  to  an  implicit  representation  in  a  regularly 
spaced  3D  grid.  The  result  is  visualized  again  as  a  triangulated  surface.  Figure  7  shows  a 
rendering  of  a  scanned  version  of  Michelangelo’s  David  which  has  several  holes. 

Figures  8(a)  and  8(b)  show  some  particular  holes  with  a  bounding  box  isolating  them. 
Figures  8(c)  and  8(d)  show  the  triangulated  surface  (the  data)  around  the  hole.  The 
reconstructed  surface  by  the  JIVFGL  model  are  shown  in  Figures  8(e)  and  8(f).  The 
reconstructed  surfaces  look  natural. 

Figure  9(a)  shows  the  hole  in  David’s  left  hand.  Figures  9(b),(c),(d),(e)  and  (f)  show 
the  corresponding  results  obtained  minimizing  the  absolute  value  of  the  Laplacian,  the 
square  of  the  Laplacian,  diffusion  of  curvature,  the  AMLE,  and  Laplace  equations,  respec¬ 
tively.  All  reconstructions  look  again  natural,  but  we  observe  that  the  result  obtained 
with  curvature  diffusion  is  less  smooth.  The  differences  between  reconstructions  are  more 
clearly  seen  in  the  orthoslices  of  Figure  10. 

Figure  11(a)  shows  the  hole  in  David’s  hair.  Figures  11(b), (c),(d),(e)  and  (f)  show  the 
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(a)  An  ideal  pyramid. 


(b)  Pyramid  with  a  gap. 


(f)  Curvature  diffusion. 


(g)  AMLE. 


(h)  Laplace  equation. 


Figure  3:  Reconstruction  of  a  pyramid  with  a  gap. 
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(a)  Pyramid  and  orthoslice  10  of  (b)  Orthoslice  of  the  ideal  recon- 


z  axis. 


(c)  Reconstruction  using 

JIVFGL. 


(f)  Curvature  diffusion. 


st  ruction. 


(d)  Minimizing  the  integral  of 


(e)  Minimizing  the  integral  of 


I  A-D|. 


|AD|2. 


(g)  AMLE. 


(h)  Laplace  equation. 


Figure  4:  Orthoslices  of  the  reconstruction  of  a  pyramid  with  a  gap. 
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(a)  An  ideal  torus. 


(b)  Torus  with  a  gap. 


(c)  Reconstruction  using 


JIVFGL. 


(d)  Minimizing  the  integral  of  (e)  Minimizing  the  integral  of 
|  AD|.  \AD\2. 


(f)  AMLE. 


(g)  Laplace  equation. 


Figure  5:  Reconstruction  of  a  torus  with  a  gap. 


(a)  Slice  of  the  ideal  torus.  (b)  The  gap  of  the  torus.  (c)  Reconstruction  using 

JIVFGL. 


(d)  Minimizing  the  integral  of  (e)  Minimizing  the  integral  of 
|  AD|.  |  AD|2. 


(f)  AMLE. 


(g)  Laplace  equation. 


Figure  6:  Orthoslices  of  the  reconstruction  of  a  torus  with  a  gap. 
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Figure  7:  Scanned  version  of  Michelangelo’s  David  (data  from  the  Stanford  Michelangelo’s 
project) 

corresponding  results  obtained  minimizing  the  absolute  value  of  the  Laplacian,  the  square 
of  the  Laplacian,  diffusion  of  curvature,  the  AMLE,  and  Laplace  equations,  respectively. 
Again,  the  comparison  between  different  reconstructions  can  be  done  with  the  orthoslices 
displayed  in  Figure  12. 

7  Conclusions 

In  this  note  we  have  shown  how  to  extend  our  previous  work  on  variational  image  inpaint¬ 
ing  to  fill-in  surface  holes.  The  idea,  inspired  by  [21]  and  [7,  8],  is  to  represent  the  surface 
of  interest  by  means  of  a  function  u,  as  an  upper  level  set  [u  >  0] ,  and  minimize  an  energy 
functional  which  integrates  a  power  of  the  mean  curvature  of  the  level  sets  of  u.  Then  we 
use  a  gradient  descent  method  and,  thus,  we  run  a  system  of  coupled  geometric  partial 
differential  equations  that  permit  to  geometrically  continue  the  surface  into  the  hole.  We 
have  also  discussed  other  curvature  based  hole  surface  reconstruction  models,  one  of  them 
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(c)  A  zoomed  detail  of 
(a)  showing  the  triangulated 
surface  with  the  hole. 


(d)  A  zoomed  detail  of 
(b)  showing  the  triangulated 
surface  with  the  hole. 


(e)  The  reconstruction  of  the 
hole  in  (c). 


(f)  The  reconstruction  of  the 
hole  in  (d). 


Figure  8:  Two  different  parts  of  David’s  surface  and  their  corresponding  bounding  boxes 
Q  containing  a  hole  on  the  surface.  The  corresponding  results  obtained  using  JIVFGL. 
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(a)  The  hole  in  David’s  left 
hand. 


(b)  Minimizing  the  absolute 
value  of  the  Laplacian. 


(c)  Minimizing  the  square  of  the  (d)  Curvature  diffusion. 

Laplacian. 


(e)  AMLE. 


(f)  Laplace  equation. 


Figure  9:  The  hole  in  David’s  left  hand  and  some  reconstructions  using  different  methods. 
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(a)  Hand  surface  and  coronal  or-  (b)  Orthoslice  of  the  gap. 

thoslice  number  35. 


(c)  Reconstruction  using  (d)  Minimizing  the  absolute  (e)  Minimizing  the  square 

JIVFGL.  value  of  the  Laplacian.  of  the  Laplacian. 


(f)  Curvature  diffusion.  (g)  AMLE.  (h)  Laplace  equation. 

Figure  10:  Orthoslices  of  the  reconstruction  of  the  hole  in  David’s  left  hand. 
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(a)  The  hole  in  David’s  hair. 


(b)  Minimizing  the  integral  of 
|AD|. 


(c)  Minimizing  the  integral  of  (d)  Curvature  diffusion. 

|  AD|2. 


(e)  AMLE. 


(f)  Laplace  equation. 


Figure  11:  The  hole  in  David’s  hair  and  some  reconstructions  using  different  methods. 
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(a)  Hair  surface  and  axial  or-  (b)  Orthoslice  of  the  gap. 

thoslice  number  59. 


(c)  Reconstruction  using  (d)  Minimizing  the  integral  of  (e)  Minimizing  the  integral  of 

JIVFGL.  |AD|.  |AD|2. 


(f)  Curvature  diffusion. 


(g)  AMLE. 


(h)  Laplace  equation. 


Figure  12:  Orthoslices  of  the  reconstruction  of  the  hole  in  David’s  hair. 
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based  on  a  variational  model  which  integrates  the  Laplacian  of  a  distance  function,  the 


other  is  heuristic  and  is  based  on  the  diffusion  of  a  function  which  represents  the  mean 
curvature  of  level  sets  of  an  underlying  implicit  function.  In  all  these  cases,  we  have  shown 
reconstruction  of  surface  holes  both  for  synthetic  and  real  data. 

Finally,  we  have  also  shown  simpler  methods  based  on  the  Laplace  equation  and  the 
so-called  AMLE  model  which  reconstructs  a  function  which  is  distance-like  near  the  known 
part  of  the  surface  and  whose  zero  level  set  can  be  interpreted  as  the  reconstructed  surface. 
If  our  interest  is  just  to  find  a  smooth  reconstruction,  this  approach  may  be  sufficient.  If 
one  wants  a  reconstruction  which  is  based  on  minimizing  mean  curvature,  it  can  just  serve 
as  an  initialization  stage. 


Acknowledgments 

This  paper  is  dedicated  to  Andres  Francisco  Sole  Martinez,  a  friend  and  colleague  who 
departed  this  world  too  early. 

GS  wants  to  thank  Prof.  Marc  Levoy  from  Stanford  University  for  providing  him  an 
early  version  of  his  paper  on  hole  filling.  This  paper  and  follow  up  conversations  with 
Prof.  Levoy  resulted  in  challenging  ourself  to  extend  our  techniques  for  image  inpainting 
to  the  problem  of  hole  filling,  leading  to  the  work  here  presented.  We  also  thank  the 
Stanford  Michelangelo  Project  for  data  provided  for  this  work,  Marcelo  Bertalnuo  who 
constantly  provides  good  feedback  to  us  on  inpainting  problems,  and  Colorna  Ballester 
for  her  help  at  several  instances  of  this  work.  The  second  author  was  partially  supported 
by  a  grant  from  the  Spanish  Ministerio  de  Ciencia  y  Tecnologfa,  reference  FP2000-5801 
and  by  the  Institute  for  Mathematics  and  its  Applications  (University  of  Minnesota). 
The  first,  second  and  fourth  authors  acknowledge  partial  support  by  the  Departament 


38 


d’Universitats,  Recerca  i  Societat  de  la  Informacio  de  la  Generalitat  de  Catalunya  and  by 


PNPGC  project,  reference  BFM2003-02125.  The  third  author  was  partially  supported  by 

the  Office  of  Naval  Research,  the  National  Science  Foundation,  the  National  Geospatial- 

Intelligence  Agency,  and  DARPA. 

References 

[1]  L.  Ambrosio,  N.  Fusco  and  D.  Pallara,  Functions  of  Bounded  Variation  and  Free 
Discontinuity  Problems,  Oxford  Mathematical  Monographs,  2000. 

[2]  N.  Amenta,  M.  Bern,  and  M.  Kamvysselis,  A  new  Voronoi-based  surface  reconstruction 
algorithm,  Proc.  SIGGRAPH  1998,  ACM,  (1998)  415-421. 

[3]  Arnira,  Arnira  Visualization  and  Modeling  System,  http://www.AmiraVis.com 

[4]  G.  Aronsson,  Extension  of  functions  satisfying  Lipschitz  conditions,  Ark.  for  Math.,  6 
(1967)  551-561. 

[5]  G.  Aronsson,  On  the  partial  differential  equation  uxuxx  +  2uxuyuxy  +  uyuyy  =  0,  Ark. 
for  Math.,  7  (1968)  395-425. 

[6]  C.L.  Bajaj,  F.  Bernardini,  and  G.  Xu,  Automatic  reconstruction  of  surfaces  and  scalar 
fields  from  3D  scans,  Proc.  ACM  SIGGRAPH  1995,  ACM,  (1995)  109-118. 

[7]  C.  Ballester,  M.  Bertalmio,  V.  Caselles,  G.  Sapiro,  and  J.  Verdera,  Filling-in  by  joint 
interpolation  of  vector  fields  and  grey  levels,  IEEE  Trans.  Image  Processing  10  (2001) 
1200-1211. 

[8]  C.  Ballester, V.  Caselles,  and  J.  Verdera,  Dissoclusion  by  joint  interpolation  of  vector 
fields  and  gray  levels,  Multiscale  Model.  Simul.  2  (2003)  80-123. 


39 


[9]  T.  Battacharya,  E.  DiBenedetto  and  J.  Manfredi,  Limits  as  p  — ►  oo  of  A pup  =  /  and 
related  extremal  problems,  Rendiconti  Sem.  Mat.  Fascicolo  Speciale  NonLinear  PDEs, 
Univ.  di  Torino  (1989)  15-68. 

[10]  F.  Bernardini,  J.  Mittleman,  H.  Rushmeier,  C.  Silva,  and  G.  Taubin,  The  ball-pivoting 
algorithm  for  surface  reconstruction,  IEEE  Transactions  on  Visualization  and  Com¬ 
puter  Graphics  5  (1999)  349-359. 

[11]  M.  Bertalmio,  A.  L.  Bertozzi,  and  G.  Sapiro,  Navier-Stokes,  fluid  dynamics,  and 
image  and  video  inpainting,  Proc.  IEEE  Computer  Vision  and  Pattern  Recognition 
(CVPR),  Hawaii,  (2001)  355-362. 

[12]  M.  Bertalmio,  G.  Sapiro,  V.  Caselles,  and  C.  Ballester,  Image  inpainting,  Proc.  SIG- 
GRAPH  2000,  ACM,  (2000)  417-424. 

[13]  M.  Bertalnuo,  L.  Vese,  G.  Sapiro,  and  S.  Osher,  Simultaneous  texture  and  structure 
image  inpainting,  IEEE  Transactions  on  Image  Processing  12  (2003)  882-889. 

[14]  J.C.  Carr,  R.K.  Beatson,  J.B.  Cherrie,  T.J.  Mitchell,  W.R.  Fright,  B.C.  McCallum, 
and  T.R.  Evans,  Reconstruction  and  representation  of  3D  objects  with  radial  basis 
functions,  Proc.  SIGGRAPH  2001,  ACM,  (2001)  67-76. 

[15]  V.  Caselles,  J.M.  Morel  and  C.  Sbert,  An  axiomatic  approach  to  image  interpolation, 
IEEE  Transactions  on  Image  Processing  7  (1998)  376-386. 

[16]  T.  Chan  and  J.  Shen,  Local  inpainting  models  and  TV  inpainting,  SIAM  J.  Appl. 
Math.  62  (2001)  1019-1043. 


40 


[17]  Y.  Chen  and  G.  Medioni,  Description  of  complex  objects  from  multiple  range  im¬ 
ages  using  an  inflating  balloon  model,  Computer  Vision  and  Image  Understanding,  61 
(1995)  325-334. 

[18]  U.  Clarenz,  M.  Rurnpf,  and  A.  Telea,  Finite  elements  on  point  based  surfaces,  Com¬ 
puters  and  Graphics,  to  appear,  2005. 

[19]  U.  Clarenz,  U.  Diewald,  G.  Dziuk,  M.  Rurnpf,  and  R.  Rusu,  A  finite  element  method 
for  surface  restoration  with  smooth  boundary  conditions,  Computer  Aided  Geometric 
Design  21  (2004)  427-445. 

[20]  B.  Curless  and  M.  Levoy,  A  volumetric  method  for  building  complex  models  from 
range  images,  Proc.  SIGGRAPH  1996,  ACM,  5  (1996)  303-312. 

[21]  J.  Davis,  S.  Marschner,  M.  Garr,  and  M.  Levoy,  Filling  holes  in  complex  surfaces 
using  volumetric  diffusion,  First  International  Symposium  on  3D  Data  Processing, 
Visualization,  and  Transmission  (2002). 

[22]  J.  S.  De  Bonet,  Multiresolution  sampling  procedure  for  analysis  and  synthesis  of 
texture  images,  Proc.  SIGGRAPH  1997,  ACM,  (1997)  361-368. 

[23]  H.  Dinh,  G.  Turk,  and  G.  Slabaugh,  Reconstructing  surfaces  using  anisotropic  basis 
functions  ,  Proceedings  IEEE  Int.  Conference  on  Computer  Vision  (2001)  606-613. 

[24]  H.  Edelsbrunner  and  E.P.  Mucke,  Three-dimensional  alpha  shapes,  ACM  Transac¬ 
tions  on  Graphics  13  (1994)  43-72. 

[25]  A.  A.  Efros  and  T.  K.  Leung,  Texture  synthesis  by  non-parametric  sampling,  IEEE 
International  Conference  on  Computer  Vision,  Corfu,  Greece,  (1999)  1033-1038. 


41 


[26]  D.  Heeger  and  J.  Bergen,  Pyramid  based  texture  analysis/synthesis,  Proc.  SIG- 
GRAPH  1995,  ACM,  (1995)  229-238. 

[27]  H.  Hoppe,  T.  DeRose,  T.  Duchamp,  J.  McDonald,  and  W.  Stuetzle,  Surface  recon¬ 
struction  from  unorganized  points,  Proc.  SIGGRAPH  1992,  ACM,  (1992)  71-78. 

[28]  K.  Ito  and  K.  Kunisch,  An  active  set  strategy  based  on  the  augmented  Lagrangian 
formulation  for  image  restoration,  M2AN  Math.  Model.  Nurner.  Anal.  33  (1999)  1-21. 

[29]  R.  Jensen,  Uniqueness  of  Lipschitz  extensions:  Minimizing  the  sup  norm  of  the  gra¬ 
dient,  Arch.  Rat.  Mech.  Anal.  123  (1993)  51-74. 

[30]  P.  Juutinen,  Absolutely  minimizing  Lipschitz  extensions  on  a  metric  space,  Annales 
Acadademiae  Scientiarum  Fennicae  Mathematica  27  (2002)  57-67. 

[31]  M.  Levoy,  K.  Pulli,  B.  Curless,  S.  Rusinkiewicz,  D.  Roller,  L.  Pereira,  M.  Ginzton, 
S.  Anderson,  J.  Davis,  J.  Ginsberg,  J.  Shade,  and  D.  Fulk,  The  Digital  Michelangelo 
Project:  3D  scanning  of  large  statues,  Proc.  SIGGRAPH  2000,  ACM,  (2000)  269-276. 

[32]  S.  Masnou  and  J.  Morel,  Level-lines  based  disocclusion,  IEEE  Int.  Conf.  Image  Pro¬ 
cessing  (1998)  259-263. 

[33]  M.  Nitzberg,  D.  Mumford,  and  T.  Shiota,  Filtering,  Segmentation,  and  Depth, 
Springer- Verlag,  Berlin,  1993. 

[34]  J.  M.  Ogden,  E.  H.  Adelson,  J.  R.  Bergen,  and  P.  J.  Burt,  Pyramid-based  computer 
graphics,  RCA  Engineer  30  (1985)  4-15. 

[35]  S.  Osher  and  R.  Fedkiw,  Level  set  methods  and  dynamic  implicit  surfaces,  Applied 
Mathematical  Sciences,  153,  Springer- Verlag,  New  York,  2003. 


42 


[36]  S.  Osher  and  J.A.  Sethian,  Fronts  propagating  with  curvature  dependent  speed:  Al¬ 
gorithms  based  on  Hamilton- Jacobi  formulation,  J.  Computational  Physics  79  (1988) 
12-49. 

[37]  S.  Rane,  G.  Sapiro,  and  M.  Bertalmio,  Structure  and  texture  filling-in  of  missing 
image  blocks  in  wireless  transmission  and  compression  applications,  IEEE  Trans.  Image 
Processing  12  (2003)  269-303. 

[38]  E.  Simoncelli  and  J.  Portilla,  Texture  characterization  via  joint  statistics  of  wavelet 
coefficient  magnitudes,  Proc.  5th  IEEE  Int.  Conf.  on  Image  Processing  (1998)  62-66. 

[39]  G.  Turk,  and  M.  Levoy,  Zippered  polygon  meshes  from  range  images,  Proc.  SIG- 
GRAPH  1994,  ACM,  (1994)  311-318. 

[40]  J.  Verdera,  V.  Caselles,  M.  Bertalnuo,  and  G.  Sapiro,  Inpainting  surface  holes,  IEEE 
International  Conference  on  Image  Processing,  ICIP  2003,  Barcelona,  Spain,  (2003) 
903-906. 

[41]  M.D.  Wheeler,  Y.  Sato,  K.  Ikeuchi,  Consensus  surfaces  for  modeling  3D  objects  from 
multiple  range  images,  Proc.  IEEE  Int.  Conference  on  Computer  Vision,  (1998)  917- 
924. 

[42]  R.  Whitaker,  A  level-set  approach  to  3D  reconstruction  from  range  data,  International 
Journal  of  Computer  Vision  29  (1998)  203-231. 

[43]  Y.  Yu,  K.  Zhou,  D.  Xu,  X.  Shi,  H.  Bao,  B.  Guo  and  H.-Y.  Shurn,  Mesh  editing  with 
Poisson-based  gradient  field  manipulation,  ACM  Transactions  on  Graphics  23  (2004) 
641-648. 


43 


[44]  H.K.  Zhao,  S.  Osher,  and  R.  Fedkiw,  Fast  Surface  Reconstruction  using  the  Level 
Set  Method,  Proc.  First  IEEE  Workshop  on  Variational  and  Level  Set  Methods,  in 
conjunction  with  Proc.  IEEE  ICCV  (2001)  194-202. 


44 


